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We construct numerical basis function sets on a lattice, whose spatial extension is scalable from 
single lattice sites to the continuum limit. 

They allow us to compute small and large bound states with comparable, moderate effort. Adopt- 
ing concepts of discrete variable representations, a diagonal form of the potential term is achieved 
through a unitary transformation to Gaussian quadrature points. Thereby the computational effort 
in three dimensions scales as the fourth instead of the sixth power of the number of basis functions 
along each axis, such that it is reduced by two orders of magnitude in realistic examples. As an 
improvement over standard discrete variable representations, our construction preserves the vari- 
ational principle. It allows for the calculation of binding energies, wave functions, and excitation 
spectra. We use this technique to study central-cell corrections for excitons beyond the continuum 
approximation. A discussion of the mass and spectrum of the yellow exciton series in the cuprous 
oxide, which does not follow the hydrogenic Rydberg series of Mott-Wannier excitons, is given on 
the basis of a simple lattice model. 



I. INTRODUCTION 

The properties of elementary excitations in solids can 
often be understood in a continuum approximation. A 
prominent example are bound electron-hole pairs in semi- 
conductors, the cxcitonsi. Mott-Wannier— excitons 
with large radius, found in materials such as silicon or 
gallium arsenide, are in first but rather good approxima- 
tion described in terms of hydrogen atoms. 

For smaller exciton radius, comparable with the lattice 
constant, deviations from the hydrogen-like properties 
occur. Central-cell corrections become important, which 
account for the possibility of finding electron and hole 
in the same unit cell and include non-parabolic disper- 
sions and modifications of the 1 / r-Coulomb potential^. A 
prototypical material is the cuprous oxide CU2O, which 
receives constant attention in the search for an excitonic 
Bose-Einstein condensato^i^. It features an interesting 
property: The mass of the "yellow" Is excitons (2.6mo, 
where toq is the free electron mass) is larger than the 
sum of the electron (LOtoq) and hole mass (0.7mo)^. The 
mass enhancement, after all by 50%, indicates the break- 
down of the continuum approximation. 

It is our intention to study central-cell corrections for 
excitons starting from a microscopic lattice model. In the 
present paper we restrict ourselves to a simple two-band 
model. Despite its limitations the discussion will provide 
us with a first understanding and set the reference for 
future work. Improved studies, e.g. with the inclusion 
of realistic ab-initio band structures, and extensions to 
biexcitons and exciton-exciton scattering are under cur- 
rent investigation. 

The study of excitons provides motivation for the de- 
velopment of numerical methods for few-particle systems 
on lattices. For small exciton radius, eigenstates can be 
obtained from a "plain" lattice calculation, where the nu- 
merical wave function is restricted to a finite number of 



lattice sites. With maximal distance L between electron 
and hole, wave function values at {2L + 1)^ lattice sites 
must be stored in memory for a three-dimensional (3D) 
problem. As a consequence of the scaling the resource 
consumption of such calculations increases rapidly with 
the exciton radius. Already for a radius of, say, 10 lattice 
sites we must deal with about 10^ wave function values. 
Excited states are obtained with even higher effort. For 
biexcitons with four particles, the scaling is (x L^. 

The above numbers indicate that a different approach 
is needed. In the present work we introduce a varia- 
tional lattice formulation of discrete variable representa- 
tions (DVR) developed for molecular physics and theo- 
retical chemistry applications^rJi. The wave function is 
represented in a product basis of sine functions, whose 
spatial width is a free parameter that is varied in order 
to minimize the energy and thus optimize the numerical 
wave function. For small basis function width, this ap- 
proach reduces to a plain lattice calculation. Allowing 
the basis function width to grow the transition to the 
continuum limit is addressed with constant effort which 
is independent of the wave function radius. 

For a decisive reduction of the computational effort 
we rely on the DVR idea of using Gaussian quadrature 
for the potential term in the Hamilton operator. With 
N basis functions along each coordinate axis in 3D, a 
straightforward variational calculation in the sine basis 
requires all N^^^ potential matrix elements between the 
A'tot = basis functions. These matrix elements have 
to be calculated and then used in each application of the 
Hamilton operator. With the DVR, the potential term is 
represented by a diagonal matrix and thus requires only 
A^tot matrix elements. Since the kinetic energy is sepa- 
rable in appropriate coordinates, the total effort scales 
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For N = 100, 



only as A^"^ = N^^^ instead of 
the effort is reduced by a factor of the order 10000. 
We deviate in two points from the standard DVR for- 
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mulation. First, we use a lattice function basis. Second, 
we adapt the DVR construction to circumvent its sin- 
gle drawback, the violation of the variational principle 
through the Gaussian approximation of the potential. In 
our formulation, the DVR strictly obeys the variational 
principle. This is important because it allows for the de- 
termination of the optimal basis function width through 
energy minimization. Since we do not know the wave 
function size, e.g. the excitonic lattice Bohr radius, in 
advance, the selection of a suitable basis function width 
would otherwise be difficult and error-prone. With these 
two modifications, our approach covers the entire range 
of small bound states occupying few lattice sites through 
the intermediate regime up to the limit of large weakly 
bound states in the continuum. 

The paper is organized as follows. In Sec. |lT]we intro- 
duce the exciton lattice model for the present investiga- 
tion, and derive the central relation between exciton mass 
and kinetic energy. In Sec. IIIII the variational sine basis 
is defined, the variational DVR is explained and applied 
to two test examples. In Sec. IIVI we study central-cell 
corrections for small radius excitons using results from 
variational DVR calculations for the lattice model. After 
a discussion of the yellow exciton series in the cuprous 
oxide, we summarize our findings in Sec. [V] 



II. THE EXCITON MODEL 

We study 3D excitons within a simple two-band model 
on a cubic lattice, with a cosine dispersion 
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for the conduction {Ef.) and valence [Eh) band. Instead 
of the electron mass me and hole mass m^ , we will also 
use the reduced (m^) and total (M) electron- hole mass 



m,. = m„ ^ 



, 7\/ = nip 
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of the kinetic energy of the electron 



Tp = 
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the kinetic energy of the hole 
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and the attractive Coulomb interaction 
U = Y,U{r-r')clc^hl,h, 
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between both, with 
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In these expressions, ci^^ and hi^^ denote fermionic op- 
erators for an electron or hole at lattice site r. The lat- 
tice sites are indexed by integer numbers, i.e. G Z 
(this explains the appearance of the parameter a in the 
Coulomb interaction) . The unit vector along each axis is 
denoted by e^. The Hamilton operator in the form given 
has five parameters, some of which are redundant as will 
be seen later: The electron/hole masses m^./f^, the (ef- 
fective) lattice constant a, the dielectric constant e, and 
the local Coulomb factor V. While e characterizes the 
long-range part of Coulomb interaction, the parameter 
V describes the relative strength of electron and hole in- 
teraction in the same unit cell. It should be ^ > 1 since 
|C/(0)|>|C/(r^O)|. 



B. Separation of center-of-mass motion 

The exciton wave function can be written as 

\i')=Y.^ir,r')clhUv^c) (8) 



The parameter a, with physical dimension 'length', deter- 
mines the typical exciton radius below which central-cell 
corrections are significant. It could be identified as the 
lattice constant for realistic band dispersion, but should 
be considered as an effective model parameter in the 
present simple treatment, similar to Ref. 0- 



A. The lattice model 

The lattice Hamilton operator is given as the sum 

H^Te + Th + U (3) 



where |vac) denotes the semiconductor ground state with 
filled valence and empty conduction band. We are inter- 
ested in wave functions with definite exciton momentum 
UK.. Translational invariance requires 



V'(r + R, r' + R) = e'"'^-^7/'(r, r') 



(9) 



which allows for separation of the center-of-mass motion 
through the ansatz 



■0(r,r') = e 
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Every combination kg + k/i = K is allowed, and any two 
choices are related by a unitary transformation. 
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The wave function of relative electron-hole motion (j){r) 
obeys the effective one-particle Schrodinger equation 



i?0(r) 



3^2 
— — < 



ia(5ke-ei „ — ia(5kh-ei 



2a^ '■ — ' \ vrif, rufi 



C/(r)<^(r) 



(11) 



In general, this equation is complex. For the simple co- 
sine dispersion assumed here, however, the choice 
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sin afc, 
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cos aki + ^ ' 



leads to a Schrodinger equation Ecj) = Hcj) with the real 
Hamilton operator 



H(l){r) 
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For small |K| ^ tt/o, we can expand the cosine and 
square root in H to second order, to find the Hamilton 
operator for low momentum states 
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(15) 

where, similar to Eqs. (HI), (IS|), the kinetic energy opera- 
tors Ti(j){r) = 2<p{r) — (j){r + e.;) — (j){r — e^) are used. 

For K = 0, the total mass M drops out of the Hamil- 
ton operator iJ<g. In this respect the lattice problem 
resembles the continuum problem, in so far as only the 
reduced mass determines the K = eigenstates and 
energies. Differences may however arise at finite K. 

We can now perform the continuum limit for by 
letting the effective lattice constant a — > for fixed values 
of the other parameters. In this limit, the wave function 
goes over into a continuous function 0c (or) = 4>{r). The 
kinetic energy operator reduces to the the derivative op- 



erator Til or 



-d /dr~, and V drops out. We obtain 



the continuum Hamilton operator 
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the Hamilton operator of a hydrogen-like atom. The 
eigenenergies for K = are En = —Rx/n^ with the ex- 
citonic Rydberg Rx, and the lowest state wave function 
radius is the excitonic Bohr radius as- Both quantities 
are given by 



2h^e^ ' 



(17) 



We choose Rx as the unit of energy and qb as the unit 
of length for the lattice problem. The exciton wave func- 
tions and energies E/Rx at K = are the eigenfunctions 
and eigenenergies of the dimensionless Hamilton operator 
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where 



w(0) = y, M(r) = l/|r| for r 7^ 



(18) 



(19) 



Only two dimensionless parameters, as /a and V, oc- 
cur. The parameter as /a distinguishes between the lat- 
tice regime as /a < 1 with significant central-cell cor- 
rections and the continuum - or Mott-Wannier - regime 
as/ a ^ 1, where central-cell corrections are absent. In 
the limit /a — > oo the eigenenergies of this Hamil- 
ton operator approach the normalized hydrogen values 
En ~ — I/ti^, independent of V. 



C. Central-cell corrections and the exciton mass 



The exciton mass is defined as 
92 



Mx = 

X QKj 



i?(K)|K=0 , 
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where E{K.) denotes the exciton binding energy as a func- 
tion of K (because of isotropy of the Hamilton operator 
the result is independent of i — x,y,z). From Eq. (|15p it 
follows that 



M 



M 



X 



l-i< 



(21) 



where (/jq (r) is the lowest state wave function of Hx ■ Note 
that for a (low-energy) bound state < (T^) < 1. In gen- 
eral, Mx/M > 1, and the exciton is heavier than electron 
and hole combined. In the continuum limit, (Tj) — s- as 
the wave function becomes larger, and Mx = M. 

Interestingly, by Eq. (|2ip the mass enhancement 
Mx I ^ depends only on the (normalized) kinetic energy 
in the lowest exciton state at K = 0. While this state- 
ment is trivially true in the continuum limit, it is some- 
what surprising for the lattice case. As discussed above, 
the K = wave function depends on as /a and V only. 
We see again that out of the five basic parameters of the 
model only two combinations determine the importance 
of central-cell corrections. 
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Based on the simple cosine dispersion, Eq. (|21|) is too 
crude to give accurate results for actual materials, but 
it can provide us with reasonable estimates. Consider 
for simplicity a 3D wave function given as a product of 
exponentials (j){x) oc exp{—a\x\/aB) along each axis (the 
estimates are independent of dimension). Then, {Tx) = 
2 — 2/ cosh(a/a_B), or Mx/M = cosh(a/aB). For as ~ 
a, we have Mx/M k, 1.5, but already for as = 4a it 
is Mx/M K, 1.03 close to one. Basically we see that 
central-cell corrections on the exciton mass are important 
for rather small, yet mobile, excitons. We expect that 
their properties depend strongly on the local Coulomb 
interaction V in Eq. ([7]). Furthermore, the stronger the 
binding the larger the mass enhancement. 

Eq. pip holds for any form of the interaction poten- 
tial U{y), but it fails in the most general case of different 
electron and hole dispersion. It has been derived as early 
as in Ref. [12] ~ where it was noted to explain the immo- 
bility of Frenkel excitons with high binding energy but 
we feel that the explanation given here is simpler. 

III. METHOD 

In a straightforward computational approach numeri- 
cal eigenfunctions of a lattice Hamilton operator such as 
in Eq. (fT4)) are restricted to a cubic box [l,i]'^ of finite 
extension. This requires storage of wave function values 
at lattice sites. The Hamilton operator is given by a 
sparse matrix, which allows for fast computations with 
iterative diagonalizationH or spectral algorithmsii. 

As discussed in the introduction, the problem with 
such plain lattice calculations is the growth of the nu- 
merical Hilbert space with i^, which prevents calcula- 
tions for weakly bound or excited states with large radius. 
Approaching the continuum limit the effort diverges al- 
though the wave function converges to a well-behaved 
continuous function. On the other hand, calculations in 
the continuum limit require some wave function repre- 
sentation at the beginning, since the natural space dis- 
cretization through the lattice is missing. 

We address these issues with a variational basis of sine 
functions. On the lattice, the ID basis functions are 

= I V^^-(^^) ifl<x<L, ^^^^ 
[ otherwise . 

The index n runs from 1 to the maximal value L. See 
Fig. [T] for a graphical presentation. From the ID func- 
tions, 3D basis functions are obtained as tensor products 
$;,„„(r) = $i(a;)$m(y)$„(z). Note that for notational 
convenience, we count x from 1 to L. In the numeri- 
cal calculations, where the potential is strongest at the 
origin, we work on the cube [—L^Vf'. 

The $„(a;) are orthonormal, with scalar product 

L 

($™, $„) ^ $™(a;)$„(x) = (5™„ , (23) 

x=l 





1 2 3 4 6 : 







'■'^O 5 10 15 20 5 10 15 20 



X X 

FIG. 1. First few basis and cardinal functions for the sine 
basis Eq. (|22|l for L — 21. Left panel: Basis functions 
$i(a;), . . . , $5(x) for L — 21. The positions of the lattice 
sites are indicated through the circles for the "l>2(a:)-curve. 
Right panel: Cardinal functions Xii^)^ • ■ • i Xsi^) for M = 6 
sampling points. 



but they are not complete since they vanish outside of 
the interval [l,i]. Bound state wave functions, which 
decay for large |a;| typically exponentially, can be repre- 
sented accurately for sufficiently large L. In our varia- 
tional calculations we allow the parameter L to grow for 
minimization of the approximation error. 
The ^n{x) have well-defined parity 

<i>„(L + l-a;) = (-l)"+i<I>„(x) . (24) 

The matrix elements of the kinetic energy operator are 
readily calculated as 

($^,r$„) = (2-2cos-^jJ^„ . (25) 

We see that the sine function are the eigenfunctions of 
the lattice 'particle in a box' problem, although this is of 
no further relevance here. Note that the scaling of the 
kinetic energy is 1/L'^, instead of 1/L for other functions 
in a box, as a consequence of the Dirichlet-like boundary 
conditions ^^^(O) = $„(i + 1) = 0. 

The continuum version of the sine basis is obtained 
from the lattice construction in the limit L — oo, if 
the coordinate x is scaled accordingly. We fix an in- 
terval [0,Lc], and write the basis functions as ^ni^) = 

L/Lc^{xL/Lc). For L — > cx), the basis functions are 
given by the sine functions 

i\f^sHn^) ifO<a;<Le, 

^nix) = <^ V ^ ' _ _ c , ^26) 

[ otherwise , 

where x now is a continuous variable G M. The orthonor- 
mality condition reads Jp " $,„(x)<I>„(a;) = Smn, a-nd the 

kinetic energy follows as —d^x^nix) = (^jr^ ^n{x)- 

We favor the sine basis over alternative choices because 
of the simplicity of all relevant expressions, the possibility 
of using fast Fourier transforms, and their equivalence to 
the established sine-DVRi^ in the continuum limit. 
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A. Discrete variable representation 

In a variational calculation, the ID wave function is 
given as a linear combination 

N 

==^0„$„(a;) (27) 

n=l 

of a finite number N of the basis functions (or in 
3D). Note that always N < L. In the limit N ^ L the 
variational basis is complete in the box, and the use of 
the sine basis is equivalent to a plain lattice calculation. 
Computational savings are expected for N <^ L. 

While the kinetic energy is diagonal in the sine func- 
tion basis, evaluation of the potential term requires mul- 
tiplication of the coefficients ...,(/>„} with a dense 
matrix V„m = ($„(a;), F(a;)$„(a;)), 

N 

(f>n = ^ Km0m ■ (28) 
rn—1 

Formally, this equation defines the projection V of V{x) 
on the variational Hilbert space spanned by the functions 
$i(a;),...,$jv(x), with 

N 

Vcf>{x) = ^{x) = 0n$„(x) . (29) 

n=l 

Note that V, in contrast to V{x), is not a multiplication 
operator in the x-eigenbasis. 

In Eq. (|28p . there are N'^ matrix elements in ID, but 
in 3D the effort grows as N^. The principal idea of the 
DVR to prevent this rapid growth is the approximate 
evaluation of the potential term at the sampling points 
of a Gaussian quadrature rule^. This is possible if the ba- 
sis functions are, essentially, polynomials in the position 
operator x. For the sine basis, we have 

$„(a;) =P„_i(£)$i(a:) , (30) 

with the transformed lattice position 

TTX 

X = COS 7 , (31) 

L + 1 ^ ' 

and the Chebyshev polynomials of second kind 

p^^,)^ Min + l)^rccosx] ^^^^ ^ _ ^^^^ 
sm[arccos x\ 

Recall that these polynomials satisfy a three-term recur- 
rence 

Po(.t) = 1,Pi(x)=2x, ^gg^ 

Pn+l{x) = 2xPn{x) - Pn~l{x) . 

The orthonormality of the basis functions Eq. pro- 
vides the discrete orthogonality relation 

5rn7i = {<S>rn{x),^n{x)) = (P,„ (5:)$1 , P„ (x) $1 ) 

2 • 2 Tl"^ n / '^^ \-n / TTX . 

= / sm P„,(cos )P„ cos ) 

L + L + 1 ^ L + r ^ L + l' 

(34) 



for the Chebyshev polynomials. 

The sampling points of Gaussian quadrature for the 
Chebyshev polynomials are the M roots 

Xk = cos , 1 < /c < M , (35) 

of the polynomial Pm{x). The roots are distinct, and 
— 1 < Xfe < 1. Note the slight technical complication 
that is given as a transformed position according to 
Eqs. (pO)) . (|3T]) . In original lattice coordinates we have 
Xfc = fc(L + l)/(M+l). 

Since the variational wave function (j>{x) in Eq. ([27|) 
is a linear combination of the first N basis functions, 
its construction through Eq. pop involves polynomials 
of maximal degree iV — 1. The wave function is thus 
uniquely specified through the values 

N 

Cfc = ^ KPn-liXk) (36) 
n=l 

at M sampling points x/. for every M > N. Note that in 
general x^ ^ Z, such that (t>{xk) itself is not defined. 

The DVR assumes that instead of multiplication with 
the dense matrix Vmn as in Eq. (|28|) the potential term 
can also be evaluated through the simpler multiplication 
of the wave function values with the potential values 
V{£,k) at M = N sampling points, i.e. 

N 

ik = J2 ~^nPn-l{xk) ~ V{xk)^k ■ (37) 
n=l 

Recall that </)„ are the coefficients of the wave function 
V4>{x) in the sine basis from Eq. (|28p . 

Eq. (|37p is exact in the full Hilbert space (i.e. for 
M — N — L), where it just states that the potential acts 
as a multiplication operator in the position eigenbasis, i.e. 
{y(j)){x) = V{x)4>{x). In general, it is an approximation 
because the projection V of the potential operator onto 
the variational Hilbert space is no longer a multiplica- 
tion operator. Since Gaussian quadrature with N points 
is exact for polynomials of maximal degree 2N — 1, the 
approximation is expected to become accurate for suffi- 
ciently large M, N . 

The benefits of the approximation Eq. ([57]) are two- 
fold: First, the potential is now given by a diagonal op- 
eration, with N'^ instead of effort. Second, instead of 
the matrix elements Vnm only the function values ^(a;^) 
are needed. The major drawback is that the approxima- 
tion violates the variational principle. 

B. Variational discrete variable representation 

The DVR uses as many sampling points as basis func- 
tions, i.e. M = N . It has been note d^^'^^ that the ac- 
curacy of the approximation Eq. (j37p can be improved 
by using larger M > N. The question then is how large 
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M has to be chosen, e.g. for a singular potential where 
Gaussian quadrature encounters difficulties. A nice re- 
sult, which seems to have been missed in the literature, 
provides a complete answer: Independent of the poten- 
tial, M = 2N — 1 suffices for the exact evaluation of the 
potential term in the DVR fashion of Eq. ((37|). 

The crucial observation is that the projected poten- 
tial V acts on the wave function in the same way as a 
polynomial in x 



Vpix) 



2N~1 

E 



VnlPn-l{x) 



(38) 



coefficients 0„ and the values at the sampling points 
are related through an orthogonal transformation with 
the matrix 



Unk = >^kPn-l{xk) 



Tmk 



M + 1^ M+1 ' 



with Ai 



Trfc 



(41) 



M+1 M + 1 



A direct calculation shows that U'^U — UU^ = 1. The 
origin of this matrix is that it contains the eigenvectors 
of X in the basis $i(a;), . . . , ^m{x) as the columns. It is 



of maximal degree 2N ^ 2, with coefficients given by 



N 



Vni = ($„(x),$i(x)^(a;)) 



2 , . Trnx , TTX . 
sm — sm — , K ) . 



(39) 



L + l 



L + 1 L + l 



This follows from comparison of matrix elements 
in the subspace basis <I'i(x), . . . , $Ar(a::). We 
first note that Km = {'^ni{x),V{x)<^n{x)) = 

{P,n{x)Pn{x)^l{x),V{x)). Thc product P,n{x)Pn{x) is 

a polynomial of maximal degree 2N — 2, which is a linear 
combination of ^0(^^)7 • ■ • ,^2Af-2(i)- It thus suffices to 
compare the 27V — 1 matrix elements {Pn{x)<^\{x) , . . .) 
for < n < 2N — 2. By definition of Vp(x), we have 

{P.a{x)^l{x),Vp{x)) = = {Pn{x)<i>l{x),V{x)), 

using the orthogonality Eq. (|34p of the P„(x) for the 
first equality. This concludes the argument. 

Since the true projected potential V can be replaced 
identically by the polynomial Vp{x) for calculations with 
the variational wave functions, the potential term can be 
evaluated exactly through Gaussian quadrature. To a 
wave function (l){x) we associate the polynomial <f>p{x) = 



p 

1-^ n. 



(x), such that (j){x) ~ (f>p{x)^i{x). Acting 



with the potential term Vp{x) now gives another polyno- 
mial Vp{x)(j)p{x) of maximal degree 3N — 3. From this, 
the polynomial (t>p{x) of the new wave function must be 
obtained, again of maximal degree — 1. In total, we 
deal with polynomials of maximal degree 4A'^ — 4. Since 
a Gaussian quadrature with M sampling points is exact 
for polynomials of maximal degree 2M — 1. the choice 
M — 2N — 1 guarantees that the expression 



^k = Vp{xk)^k 



(40) 



is exact. Wc can thus evaluate the potential term ex- 
actly by using twice as many sampling points as basis 
functions, and the values Vp{xk) instead of the potential 
function values V{xk)- This provides us with the effi- 
ciency benefits of thc DVR and preserves the variational 
principle. 



C. Implementation 

Let us now explain how we use the variational DVR 
(VDVR) in practice. First note that the wave function 



(42) 



Note that for AI > N only a rectangular N x M sub- 
matrix of Unk is used. Note further that only L ba- 
sis functions exist on the interval [1 , i] , such that M = 
min{2A- 1,L}. 

As a side remark, we mention the eigenfunctions 
Xk{x) = Yln=iUnk^n{x) of x in the Hilbert space 
spanned by ^jv/, the cardinal functions. They 

can be used to represent (j){x) = X]a;=i ^k£,kXk{x) directly 
through the values at the sampling points. Typical 
cardinal functions are depicted in Fig. [1] For thc maxi- 
mal value M = L, Xkix) = 5xk is the lattice (5-function 
localized at the single site x = k. We here recover the 
plain lattice calculation. We will not use cardinal func- 
tions explicitly in this work. 

We can now proceed as follows for thc evaluation of 
the potential term: Obtain the values S,k from the coef- 
ficients (j>n through transformation with Unk (Eq- ((42|) ). 
then multiply each with Vp{xk) according to Eq. (jiU)) . 
and transform back to the new coefficients d>„. That is, 



M 
k=l 



N 



Vpixk) Ur, 



(43) 



or more concisely (j) — UVpU^cj) where the matrix Vp = 
{Vp{xk)Sjk)jk is diagonal. Note that, in contrast to the 
standard DVR, we consider the 0„ instead of the £,k as 
the primary objects in the calculation since it simplifies 
the formulation and calculation for M > N (cf. Ref . [iTf) . 

The ID Hamilton operator in the (V)DVR formulation 
is given by 



VDVR 



= UVpU+ + T 



(44) 



with thc diagonal N x N kinetic energy matrix T from 
Eq. p5|) . thc diagonal M x AI potential matrix Vp, and 
the rectangular N x M transformation matrix U. 

The extension to 3D is straightforward. The Hamilton 
operator retains the form of Eq. . The kinetic energy 
is the sum + Ty + Tz and remains diagonal, and the 
diagonal potential matrix has now entries. The trans- 
formation matrix is a tensor product C/'-'^' = U ®U ®U , 

which acts as ^^fc = Y^lmn UtlUjmUkn(f>lmn- 
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The multiphcation with the f/^"^^ matrix is best done 
sequentially. Since in the tensor products each matrix 
applies only along a single axis, every multiplication re- 
quires operations. The matrices V and T are diag- 
onal, and require operations. The total operation 
count is thus of the order iV^ instead of N^. A certain 
overhead for the additional transformations is present. 
With a more detailed counting we find that already for 
= 10 the effort is reduced by a factor 3 in compari- 
son to a non-DVR evaluation of the potential term. For 
N = 100, the reduction is by a factor of 360: A calcula- 
tion that takes ten minutes with the (V)DVR would take 
two and a half days without! The numbers become even 
more favourably if we use a fast Fourier transform (FFT) 
for the multiplication with t/'-'^^. 

D. Calculation of Vp{xk) 

The potential enters Eq. gO]) through the M = 2A - 1 
matrix elements Vp{xk)- They are obtained from the Vni 
in Eq. p9p through transformation with Unk as 

2JV-1 ^ 2N-1 

Vp{xk) = Pn-l{S:k)Vnl ^ — UnkVnl ■ (45) 

In practical applications the evaluation of the scalar 
product in Eq. (j39p . which requires summation over L 
(or L^) lattice sites, is not desirable. If the potential is 
given by a smooth function V{x) we can use Gaussian 
quadrature instead. With Ng > M sampling points, we 
obtain the approximation 

n — 1 j — 1 

where U^^ , A^" are defined as in Eq. (|4ip with Ng replac- 
ing M. This expression essentially describes the projec- 
tion from polynomials of maximal degree Ng onto poly- 
nomials of maximal degree M, all given through their 
values at certain sampling points. Note that the argu- 
ment of V{x) is not necessarily an integer. Therefore, we 
must know V{x) for continuous x, not only at the lattice 
sites. For Ng — Af , the Gaussian approximation reduces 
to the DVR-like expression 

Vpixk) V[kp;^) , fc = 1, . . . , M . (47) 

The entries of Vp in Eq. (jH]) can be obtained through 
any of the expressions Eqs. (|45p . pB]). (jTZ]). In practice, 
we use the simplest approximation Eq. (|47p for regular 
potentials, e.g. (an-) harmonic oscillators. For singu- 
lar potentials, such as the Coulomb potential, wc use 
Eq. (P^]) with Ng = 4 . . . 8A^, avoiding sampling points 
at the potential singularity. On the lattice we treat 
localized ((5-function) contributions, e.g. from the V- 
term in Eq. ([7]), exactly through Eq. (|45p . and use again 
Eqs. (|46p . (|T7|) for the remaining long-range part of the 
potential. 



E. Discussion 

An important conceptual difference between the 
VDVR and the original DVR formulation is the sepa- 
ration of the calculation of potential matrix elements 
through Eq. (gS]) and their actual usage in Eqs. (02]), (011). 
Two sources of error exist in the DVR: First, the re- 
placement of the full, dense matrix-vector multiplica- 
tion Eq. (1211) through the diagonal expression Eq. (j571) . 
This error is completely eliminated in the VDVR 
Eqs. (gni), (031) through the choice oi M = 2N - 1 sam- 
pling points. Second, the error incurred through approx- 
imate evaluation of the matrix elements Vp(xfc). This er- 
ror can be made smaller through better approximations 
to Eq. (051), e.g. by increasing Ng in Eq. (|46l) . 

In the original DVR, both errors contribute equally. 
The obvious way for error reduction is to increase N, 
which inflicts a computational overhead on the entire cal- 
culation. In the VDVR, the effort for a better calculation 
of the Vp{xk) needs to be invested only once in Eq. (|45l) 
or Eq. (051) before the actual use of the Hamilton oper- 
ator Eq. (0^ . but there is no reason for increasing M 
beyond 2N — 1 for the evaluation of the potential term. 

In our opinion this is a central advantage of the formu- 
lation chosen here: The calculation of matrix elements of 
V{x) is completely independent of the Gaussian quadra- 
ture underlying the diagonal VDVR evaluation of the po- 
tential term. In particular for singular potentials a large 
number Ng of sampling points in a Gaussian quadra- 
ture Eq. (j46p can be necessary to obtain accurate ma- 
trix elements, while the variational wave function is a 
good approximation already for A^, M <C Ng. It can 
also be useful to calculate the Vp{xk) with other in- 
tegration/summation procedures, e.g. adaptive Gauss- 
Kronrod-integration or specialized routines for functions 
with an integrable singularity. In the original DVR, re- 
placing the V{xk) in Eq. ((571) by better matrix elements 
is not possible, and increasing N is the only possibility 
for improvement. 

F. Examples 

1. Harmonic oscillator 

The difference between the variational and non- 
variational DVR is apparent for the (continuum) har- 
monic oscillator H ~ —^d^x + 52;^, with eigenvalues 
E„ = n — 1/2 (note that we start counting with n = 1). 
In Fig. 121 wc show the numerical result for the energy 
El, . . . , E4 of the four lowest eigenstates under variation 
of the parameter Lc- Note that the Gaussian approxima- 
tion Eq. (|47p for the potential is used. The error has two 
sources: The domain truncation error because the varia- 
tional wave function vanishes outside an interval of length 
Lc; and the basis error from the approximate represen- 
tation of the wave function through a flnitc number of 
basis functions. Since the VDVR is variational by con- 
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FIG. 2. Numerical energies of the four lowest eigenstates of 
the harmonic oscillator H = —dxx/2 + a;^/2, as a function of 
the interval length Lc for the interval [—Lc,Lc]. Shown are 
results for the VDVR (left panel) and the DVR (central and 
right panel), for A'' = 10, 21 basis functions as indicated. The 
exact eigenenergies En ~ n — 1/2 are indicated by dashed 
horizontal lines (note that we start counting with n = 1). 
The Gaussian approximation Eq. (|47|l for the potential matrix 
elements has been used. 

N El E5 Eio E20 

5 0.5001685777 5.088729909 — — 

10 0.5000004143 4.510831820 11.14028090 — 

20 0.5000000000 4.500000074 9.500334454 23.96785196 

30 0.5000000000 4.500000000 9.500000005 19.51678526 

40 0.5000000000 4.500000000 9.500000000 19.50000254 

50 0.5000000000 4.500000000 9.500000000 19.50000000 

TABLE I. Convergence of the numerical energies with the 
number of basis states N in a, VDVR calculation, for the 
eigenenergies Ei, E^, Eio, E20 of the harmonic oscillator H = 
—dxxj'i + 12 as in Fig. [S] The Gaussian approximation 
Eq. (|47|) for the potential matrix elements has been used. 
Note that the energy minimum is obtained for different opti- 
mal L (not shown), in particular for small N . 



struction, the numerical energies approximate the true 
energies from above (left panel). The optimal value is 
found by minimization of the respective numerical energy 
under variation of Lc- The values in Table U demonstrate 
the attainable precision. With only = 30 the first 10 
eigenvalues are converged, and even the 20th eigenvalues 
E-iQ is significant with a relative error below 10"'^. 

The standard DVR violates the variational principle 
even in this simple example (central and right panel in 
Fig. [2]), and it is not recovered for large N . Neverthe- 
less the DVR provides meaningful results if the plateau 
region is identified, where the DVR energy is almost con- 
stant under variation of Lc and a good approximation 
to the true energies. But the violation of the variational 
principle complicates the application of the DVR for our 
purposes, since an automatic identification of the plateau 



is difficult and error-prone. Also, we normally have no 
a-priori estimate for a suitable Lci since we do not know 
the wave function in advance. 



2. Hydrogen atom 

The second example concerns already the numerical 
calculations for the 3D cxciton problem. We consider the 
dimensionless Hamilton operator Eq. (jlSp for as /a ^ 1, 
where plain lattice calculations become increasingly de- 
manding. In Table |ll] we show the energy of the lowest 
state in the continuum limit as/a = 00, i.e. for a 3D hy- 
drogen atom, obtained with VDVR. The convergence to 
the exact value is not as favourable as for the harmonic 
oscillator (cf. Table |I| , as a consequence of the singu- 
lar potential. The singularity of the potential increases 
the error of the Gaussian approximation Eq. (|47| . The 
convergence improves for more accurate potential matrix 
elements, obtained with Eq. (|46| for Ng = 8N. Since 
the additional effort for a better calculation of the ma- 
trix elements is only invested once, prior to the actual 
use of the VDVR Hamilton operator in an iterative diag- 
onalization procedure, the running time of calculations 
increases only marginally (about 5%). 

A second error source is intrinsic to the sine basis con- 
struction, which has difficulties to resolve the cusp of 
the hydrogen wave function at r = 0. Recall that we 
approach the continuum limit starting from a cubic lat- 
tice. In the continuum limit, rotational symmetry allows 
for the separation of radial and angular coordinates and 
the choice of a better basis, e.g. of Laguerre polynomials. 
Starting from the lattice, this is prevented by the reduced 
lattice symmetry. Convergence of expansions of a rota- 
tionally invariant function with a cusp in a basis without 
rotational symmetry is relatively slow. Since full rational 
symmetry is restored only in the continuum limit, there 
is no easy fix to this problem. For smaller as /a, away 
from the continuum limit, the error is reduced. 

Despite these complications, the error with N ^ 100 is 
smaller than 3 x 10""'. The variational Hilbert space has 
dimension 10^. As a consequence of parity symmetry, 
only half of the basis functions along each axis are used, 
reducing the necessary effort by a factor of 2^ = 8. Mak- 
ing use of the full lattice symmetry requires manipulation 
of only « 21000 states for N = 100 (effectively, 30 per co- 
ordinate axis). There is probably room for improvements 
of the basis, but we are not aware of a simple solution 
to the lattice/continuum symmetry mismatch problem. 
With the present construction, we achieve an error below 
lO"'^ in all quantities and for all parameters, at modest 
computational effort. For — 30 . . . 50, the dimension 
of the variational Hilbert space is 125000 at maximum 
without consideration of lattice symmetry. Such calcula- 
tions, in double real precision, require less than 25 MB of 
computer memory. Implementing the full lattice symme- 
try reduces these numbers to 2600 states and about 2 MB 
of central storage. Recall that also the running time of 
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N 




El {Ng = 8N) 


10 


-0.92931387 


-0.96164993 


20 


-0.97321483 


-0.99026343 


30 


-0.98536003 


-0.99589247 


40 


-0.99059861 


-0.99781906 


50 


-0.99338044 


-0.99871312 


60 


-0.99492841 


-0.99912287 


70 


-0.99605729 


-0.99938235 


80 


-0.99683500 


-0.99954469 


90 


-0.99739575 


-0.99965255 


100 


-0.99781451 


-0.99972716 



exact -1.0 

TABLE II. Convergence of the ground state energy Ei with 
the number of basis states N in a, VDVR calculation, for the 
Hamilton operator (jlSp in the continuum limit as /a — oo, 
i.e. a 3D hydrogen atom. The values in the 2nd column 
are obtained with the Gaussian approximation Eq. H47p . the 
values in the 3rd column from Eq. (|46p with Ng — 8N. 



a program is significantly reduced due to the favourable 
(V)DVR scaling oc N^. For the calculations reported 
here, this level of accuracy and efficiency is sufficient. 

The real benefits of the VDVR over a plain lattice cal- 
culation become apparent in Fig. (|3]), where we show 
the energy of the lowest and first excited state of Hx 
(Eq. (HH])) as a function of as/a. We compare the VDVR 
with TV = 30 basis states per direction to a plain lat- 
tice calculation on a cube [— L, L]'^ with fixed L = 30. 
For small as /a,, the wave function radius is small and 
the VDVR reduces to the plain lattice calculation. As 
the wave function radius grows with as/a, the domain 
truncation error of the plain lattice calculation becomes 
severe and renders the results meaningless. In partic- 
ular excited states are not accessible. On the other 
hand, the error of the VDVR is bounded independently 
of the actual extension of the wave function, and the 
energies converge to the correct values Ei/Rx = — 1, 
E2/RX = —0.25 for as/a — )• 00. The error in this limit 
can be deduced from Table [TTl 

It should be noted that the optimal L in the VDVR is 
different for the lowest state and excited states, reflecting 
the larger radius of the latter. It is the advantage of the 
variational procedure to adapt itself to these differences. 



IV. EXCITONS 

According to Eq. (|2T|) and the simple estimate we gave 
in Sec. |TT1 central-cell corrections of the exciton mass be- 
come important if the exciton Bohr radius is of the or- 
der of the (effective) lattice constant, i.e. for as/a < 1. 
In Fig. |4] we show the exciton binding energy —E/Rx 
and the mass enhancement Mx/M for the lowest exci- 
ton state, as obtained with the VDVR applied to the 
Hamilton operator Eq. ^TS\\ . As we discussed in Sec. HIl 
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FIG. 3. Energy of the lowest (-Bi) and first excited state 
(E2) for the Hamilton operator Eq. (fT8)l (with V = 1), as 
a function of as I a. The solid curves are calculated with the 
VDVR and A'^ = 30 basis states. The dashed curves have been 
obtained with a plain lattice calculation on a cube [— L, I/]'^ 
with L = 30, 50. 
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FIG. 4. Exciton binding energy E / Rx and mass enhance- 
ment Mx/M as a function of as /a, for different V as indi- 
cated. Results were obtained with a VDVR calculation for 
A = 20 ... 50 with an error below 10^'^ for both quantities. 

these curves depend only on the parameters V and as/a- 
As expected, —E/Rx,MxlM — >• 1 in the continuum 
limit as/a — > cxj, independent of V. Deviations arise 
for smaller as jo,. It may be interesting to note that the 
binding energy changes more pronouncedly with y, while 
the mass enhancement, as a function of as/a, remains 
similar. The numerical data show that two opposite sit- 
uations are possible: A small binding energy and large 
mass enhancement for small V and aB/a^ and a large 
binding energy and small mass enhancement for large V 
and moderate as /a- Physically, the magnitude of both 
parameters 0,3/0- is related, e.g., to the extension of 
Wannicr functions for the conduction and valence band. 
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FIG. 5. Spectrum of the exciton model Eq. (fT8)l . for V = 3, 
as/ a = 2 (left panel) and V = 5, as /a = 1 (right panel), cal- 
culated with VDVR. Grey dashed lines indicate the Rydberg 
series En/Rx = — 



For reasonable parameters a large exciton mass coincides 
with a higher binding energy, and vice versa. 

In Fig. [S] we show the typical change of the exciton 
spectrum (at K = 0) in comparison to the hydrogenic 
Rydberg series En = —Rx/n^, which is realized in the 
continuum limit as/a ^ oo. Starting from there, the n^- 
fold degeneracy of the hydrogen eigenstates is lifted in the 
lower crystal symmetry. While the notation in the figure 
refers to the hydrogen problem, a group-theoretical clas- 
sification is possible with the irreducible representation 
of the point group Oh of our lattice model, the symmetry 
group of a cubei^. The even parity "s" states (odd parity 
"p" states) arise from hydrogen states with angular mo- 
mentum I = {I = 1), and correspond to one-dimensional 
(three-dimensional) irreducible representations. For the 
even parity "d" states, the 2Z -I- 1 = 5 dimensional ir- 
reducible representation of the full rotation group splits 
into a two and a three dimensional representation under 
the reduced symmetry of Oh- Numerically we see indeed 
that each such state is a doublet of a two and three- 
fold degenerate eigenenergy, but the splitting of the order 
10^"* is not resolved in Fig. [71 

The wave functions given in Fig. [5] still resemble hy- 
drogen wave functions, although their properties, i.e. the 
binding energy or mass, do not. In our simple model, the 
energy shifts are induced by the V-tevm in Eq. ([7]), which 
affects the "s" states strongly since the probability |0(O)p 
of electron and hole being in the same unit cell is large. 
The "p,d,. . ." states are much less affected (for hydro- 
gen states, 0(0) = exactly for ^ > 1). Significant energy 
shifts thus arise from lifting of the dynamical degeneracy 
of the l/r-Coulomb potential with respect to the angular 
quantum number I, and are not associated with the split- 
ting of states with different magnetic quantum number 
m that is predicted by the lower crystal symmetry. 

To understand the significance of these results for the 




FIG. 6. Wave function (j}{x,0,0) for parameters V = 2.75, 
as/a = 0.98 as used for the cuprous oxide in Fig. [7] below. 
Shown is the lowest exciton state (left panels, with linear and 
logarithmic ordinate), the second even (top right panel) and 
the first odd state (bottom right panel). The circles give the 
values <l){x,0,0) at the respective lattice site, the solid red 
curves show the corresponding hydrogen wave functions for 
comparison. Also given is the wave function radius R. 



cuprous oxide CU2O, we show a model calculation for the 
yellow (ortho-) exciton series in this material in Fig. [Bin 
comparison to experimental data from Refs. [l^ and l20l. 
The experimental absorption energies of the odd parity 
states can be fitted to a perfect Rydberg series E^^ = 
Ef - R?^lr?, with Ef = 2.17eV for the gap energy 
and i?^ = 98.4meV for the excitonic Rydberg. With 
these values, we obtain the left panel in Fig. [7] for the 
normalized energies E j Rx=[E^^ —E^^'^) j . Obviously, 
the energies of the even parity states differ significantly 
from the Rydberg energies — 1 jv? . 

For the model calculation, we choose the two parame- 
ters V = 2.75 and as/a = 0.98 according to a fit of the 
numerical binding energy and mass enhancement of the 
lowest (Is) exciton state, as given in Fig. 01 to the exper- 
imental values E = 139meV and Mx/M = 1.5. Noth- 
ing was assumed or adjusted for higher exciton states. 
Both parameter values have a reasonable order of magni- 
tude. Recall that as /a < 1 is an indication of significant 
central-cell corrections. 

The central panel in Fig. [71 shows the calculated ex- 
citon spectrum for the above model parameters. The 
spectrum is qualitatively correct, and we find excellent 
quantitative agreement for the 2s state which reproduces 
the experimental energy with an error of only 5%. De- 
viations occur for higher even parity states, which shows 
the limitations of the simplistic model used here. 

The state labelled (IG) in the left panel does not fit 
into the model spectrum, and should probably be incor- 
porated into the green exciton series shown in the right 
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FIG. 7. Spectrum of the exciton model Eq. (|18p in compar- 
ison to experimental values for the yellow (ortho-) exciton 
series in cuprous oxid o^^'^° . The model parameters V — 2.75, 
as /a = 0.98 are determined from the binding energy and 
mass enhancement of the lowest (Is) exciton state. Shown 
are even and odd parity states for the yellow (left panel) and 
green (right panel) exciton series, and the results of the VDVR 
model calculation for the yellow series (central panel). Note 
the broken energy axis. The grey dashed lines give the ener- 
gies — of the Rydberg hydrogen series. The "IG" line in 
the left/right panel is at identical energy. 

panel in Fig. [T) Again from a fit of the experimental ener- 
gies of the odd parity states to a Rydberg series, resulting 
in = 2.31eV and = 151meV, the (IG) state is 
found exactly at the energy E^"^ — of the lowest (Is) 
green exciton state. There is further experimental ev- 
idence about the assignment of states to the yellow or 
green exciton series, such as response to strain and mag- 
netic fieldsi^i^, but with the present simple model we 
are unable to provide further analysis. 



V. CONCLUSIONS 

In the present paper we introduce a variational discrete 
variable representation for bound states on a lattice and 
apply it for a study of excitons with significant central- 
cell corrections. 

In the VDVR wave functions arc given in a variational 
basis of sine functions. It combines (i) accuracy because 
of the use of exact matrix elements and the variational 
determination of the optimal basis function width, with 
(ii) efficiency, since it evaluates the potential term in the 
DVR spirit through a diagonal matrix. We adapted the 
original DVRs in two aspects: Our construction (iii) fully 
preserves the variational principle, and (iv) bridges the 
gap between lattice and continuum calculations in a sin- 
gle unified framework. 

The example of central-cell corrections for excitons 
provides the physical motivation for the present work. 



From the simple two-band lattice model adopted here we 
can mainly draw qualitative conclusions. In the regime 
as — 0, where central-cell corrections become important, 
excitons are still closer to Mott-Wannier excitons than 
to Frenkcl excitons. Their properties, however, deviate 
significantly from the hydrogen picture - even for excited 
states where the radius exceeds the lattice constant. 

Despite the simplicity of the two-band model, wc can 
successfully reproduce the experimental spectrum of the 
yellow exciton series in the cuprous oxide CU2O, even 
with quantitative agreement. Only two parameters enter 
the model calculation, which arc fixed by the binding en- 
ergy and mass of the yellow Is exciton state. That the 
spectrum of excited states can be reproduced from two 
elementary properties of the lowest exciton state - one 
being the exciton mass with no apparent relation to ener- 
gies of optical transitions - shows how the use of a lattice 
model allows us to connect different properties through a 
fundamental microscopic description. The additional ex- 
perimental information from measurements of the exci- 
ton mass can thus be used in a theoretical interpretation 
of the exciton spectrum. 

The present model calculation is too simplistic to cover 
all relevant aspects of exciton formation. The study of 
many important effects, such as the spin-dependent en- 
ergy splitting between ortho- and para-excitons or the 
infiuence of electron-hole exchange interaction on the ex- 
citon rnass^^, has to be postponed to a forthcoming pub- 
lication, where we will also discuss how the large central- 
cell corrections for yellow excitons in comparison to the 
negligible corrections for the green exciton series are re- 
lated to the respective valence band dispersion. 

We will consider extensions of the present work in three 
directions. First, refinements of the two-band model are 
necessary. As a few principal issues, we can list (a) the 
inclusion of realistic band structures, e.g. from ab-initio 
calculations^, (b) full consideration of lattice symme- 
tries, which is particularly important for the classification 
of excited states^^, (c) improved matrix elements for the 
short-range Coulomb interaction, which can be princi- 
pally obtained from the Wannier functions of conduction 
and valence bands, and (d) the corrections to the dielec- 
tric constant for screening at short distance a^^'^^ . The 
exciton spin configuration is relevant in (e) the exchange 
interaction, leading to a splitting of exciton states^^, and 
(f) spin-orbit coupling. For comparison with experiment 
it is desirable to allow for (g) external (magnetic) fields 
and (h) strain/lattice deformatio n^ ^ 1 ^ . Also such a re- 
fined exciton model can be studied within the VDVR. 

More demanding would be the inclusion of dynami- 
cal screening, which is possible within a Green function 
formalism. We have discussed elsewhere the use of poly- 
nomial bases for Green function calculations^^, but the 
combination with VDVR is not worked out. It would give 
a polynomial basis construction both for position and en- 
ergy. Note that the present work can be understood as 
the solution of the Bethe-Salpeter equation in the special 
case of a non-frequency dependent interaction. We see no 
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urgency to include dynamical screening into the model, 
because its effect is less significant than the corrections 
listed above. 

Second, our derivation of the VDVR generalizes to ar- 
bitrary basis sets of orthogonal polynomials with only mi- 
nor modifications. In particular in the continuum limit, 
where we have more freedom for the basis choice, such a 
generalization display its full strength in comparison to 
artificial ad-hoc discretizations of position or momentum 
space. It will be discussed elsewhere. 

Third, the VDVR is powerful enough to allow for the 
study of biexcitonic systems and exciton-exciton scat- 
tering. Our discussion of the exciton spectrum shows 
why lattice models are important for an understanding 
of small-radius excitons, and the same is true for interact- 
ing two-exciton systems. A calculation of the central-cell 
corrections for exciton scattering lengths is of immediate 
relevance for Bose-Einstein-condensation. 



Bearing in mind that the cuprous oxide is one mate- 
rial where the search for Bose-Einstein-condensation of 
excitons is justified, microscopic studies of excitons in 
this material always provide us, apart from our genuine 
interest in excitons away from the Mott-Wannier limit, 
with a perspective on the conditions and limitations of 
collective exciton behaviour. In this sense, the VDVR 
technique and the lattice calculations reported here are 
one building block for the understanding of recent and 
future experiments on the cuprous oxide and similar ma- 
terials. 
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